*#第十一章 Stata编程的基础准备
*##11.1 do文档编写
*###11.1.1 do文档示例

local ef "./ef.xlsx"                                                           
import excel "`ef'", first  
foreach v of varlist _all {  
    local lab1 = `v'[1]  
    local lab2 = `v'[2]  
    local lab12 `"`lab1':`lab2'"'  
    label var `v' `"`lab12'"'  
}  
drop in 1/2  
destring, replace  

*###11.1.2 do文档编写中的Stata基本元素

search label

*局部宏(Local macros)

local pid year  
local year = 2000+21  
local n: word count learning Stata in Rstata 

/////////

display "`pid'"  //在`pid'外加上双引号，宏中存放的对象将以字符型展示
display `year'  
display `n' 

/////////

local pid year  
local year = 2000+21  
display "`pid'"  
display ``pid'' 

*全局宏(global macros)

global pid year  
global year = 2000+21  
display "$pid"  
display $year  
display $$pid  

*标量（scalar）

scalar a = 2  
scalar b = a+3  
scalar root2 = sqrt(2)  
scalar im = sqrt(-1)  
scalar s1 = "hello world"  
scalar s2 = word(s1,1)  
display a  
display b  
display root2  
display im  
display s1  
display s2  

*矩阵(matrix)

matrix input  A = (1,2\3,4) 
matrix input  B = (1,2\3,4) 
matrix define C = A + B

/////////

matrix list A   
matrix list C  

/////////

matrix A = A/A[1,1]   
matrix B = A["weight","displ"]

/////////

matrix B = I(6)   
matrix C = B[2..4, 3..6]   
matrix D = B[2..., 3...]   
matrix E = B[., 3..4] 

/////////

matrix B[2,2] = 5   
matrix list B   
matrix B[2,2] = J(3,3,5) //修改多个元素时，仅需指定修改的起始行列位置  
matrix list B  

*##11.2 编写ado命令
*###11.2.1 什么时候需要编写Stata-ado命令？

local ef "./ef.xlsx"                
import excel "`ef'", clear                                                 
labone, n(2 3) concat(:)  

*###11.2.2 ado命令的基本结构

program define helloworld   
  version 16   
  display "Hello, World."   
end  

*###11.2.3 如何编写ado命令？

quietly summarize X, detail   
local p7525 = r(p75) - r(p25)   
local p9010 = r(p90) - r(p10)   
di "p7525:" `p7525'   
di "p9010:" `p9010' 

/////////

program pctrange, rclass   
    version 14   
    syntax varlist (max=1 numeric)   
    quietly summarize `varlist', detail   
    local p7525 = r(p75) - r(p25)   
    local p9010 = r(p90) - r(p10)   
    di "p7525:" `p7525'   
    di "p9010:" `p9010'   
    return scalar p7525 = `p7525'   
    return scalar p9010 = `p9010'   
end  

/////////

sysuse auto, clear   
pctrange price   
scalar p7525=r(p7525)   
scalar p9010=r(p9010) 

/////////

cap program drop pctrange           //清除内存中可能存在的pctrange命令                                        
program pctrange, rclass   
   version 14   
   syntax varlist(max=1 numeric)[,noPRINT]   
   quietly summarize `varlist', detail   
   local p7525 = r(p75) - r(p25)   
   local p9010 = r(p90) - r(p10)   
   if "`print'"!="noprint"{   
            di "p7525:" `p7525'   
            di "p9010:" `p9010'   
    }   
    return scalar p7525 = `p7525'   
    return scalar p9010 = `p9010'   
	end  

/////////

cap program drop pctrange
program pctrange, rclass   
    version 14   
    syntax varlist(max=1 numeric) [if] [in][,noPRINT]   
    marksample touse   
    quietly summarize `varlist' if `touse', detail   
    local p7525 = r(p75) - r(p25)   
    local p9010 = r(p90) - r(p10)   
    if "`print'"!="noprint"{   
            di "p7525:" `p7525'   
            di "p9010:" `p9010'   
   }   
    return scalar p7525 = `p7525'   
    return scalar p9010 = `p9010'   
end  

/////////

pctrange price if foreign==1   
pctrange price if foreign==0

/////////

program pctrange, rclass byable(recall)   
   version 14   
   syntax varlist(max=1 numeric) [if] [in][,noPRINT]   
	marksample touse   
    quietly summarize `varlist' if `touse', detail   
    local p7525 = r(p75) - r(p25)   
	local p9010 = r(p90) - r(p10)   
    if "`print'"!="noprint"{   
            di "p7525:" `p7525'   
            di "p9010:" `p9010'   
    }   
	return scalar p7525 = `p7525'   
	return scalar p9010 = `p9010'   
	end    
	
by foreign: pctrange price 

*##11.3 mata函数编写
*###11.3.2 Mata基本语法

*Mata环境

mata:  
	 a = "Hello, World!"  
     a  
end

*Mata矩阵的相关操作

	mata:  
    A = (1,2,3 \ 4,5,6 \ 7,8,9)  
	end

/////////

	mata:  
    a=1::3  
	J(5,5,a)  
	diag(a)  
    end
	
/////////

mata:  
    D=(1,2,3 \ 4,5,6 \ 7,8,9)  
    D[1,2]  
    D[.,2]  
    D[1,.]  
    D[1::2,1]  
	D[1,1..2]  
end  

/////////

mata:  
    A = (1, 2, 3)  
    B = (4 \ 5 \ 6)  
    A+B'  
    A-B'  
    A*B  
    A:*B'  
    A:-100  
end

*Mata流程控制

mata:  
   for(i=1;i<=10;i++){  
      printf("%g squared is %g. \n",i,i^2)  
   }  
end 

/////////

mata:  
    for(i=10;i>0;i=i-2){  
        printf("%g squared is %g. \n",i,i^2)  
    }  
end  

/////////

mata:  
    i=10  
    while(i>0){  
        printf("%g squared is %g. \n",i,i^2)  
        i=i-2  
    }  
end  

/////////

mata:  
    for(i=10;i>0;i=i-2){  
        if(i==8) continue //当i==8时不执行以下循环代码，进入下一循环  
        printf("%g squared is %g. \n",i,i^2)  
    }  
end 

/////////

mata:  
    for(i=10;i>0;i=i-2){  
	        if(i==8) break //当i==8时结束循环  
        printf("%g squared is %g. \n",i,i^2)  
    }  
	end  


*###11.3.3 mata函数的基本结构

sysuse auto, clear  
putmata price = price //将price变量数据放到mata环境中的price向量
mata: sum(price)

/////////

mata:  
real vector cusum(real colvector x)  
{  //以左花括号"{"表示该函数具体代码的开始  
    n=rows(x)  
    y=J(n,1,0)  
    y[1,1]=x[1,1]  
	    for(i=2;i<=n;i++){  
	        y[i,1]=y[i-1,1]+x[i]  
	    }  
    return(y)  
	}  //以右花括号"}"表示该函数代码的结尾  
end  
mata: a = 1::10    // a是一个从1到10的列向量  
mata: b = cusum(a) // 计算a的累积加总，将结果赋予向量b  
mata: b

/////////

mata:  
void function rmean(real vector x,string scalar z){  
    n=rows(x)  
    y=0  
    for(i=1;i<=n;i++){  
        y=y+x[i]  
    }  
    y=y/n  
    st_numscalar(z,y)  
	}  
end  
	mata: a = 1::10    //a是一个从1到10的列向量  
	mata: rmean(a,"b") //计算a的均值，并将结果返回到Stata的Scalar b  
	display b  

*#第十二章 DEA模型编程的基本思路
*##12.3 步骤2：求解单个线性规划问题

	mata:   
       datatable=( 5,5,8,34,9      ///  
                \ 6,9,6,67,6       ///  
                \ 7,7,87,32,19     ///  
                \ 23,14,45,11,31   ///  
                \ 67,34,12,37,29 )  
                    
        lp = LinearProgram()               // 定义线性规划类               
        N1=5                              // 参照决策单元的个数          
        X=datatable[1..3,1]               // 取出第一个决策单元的投入数据  
	    Y=datatable[4..5,1]               // 取出第一个决策单元的产出数据  
        Xref=datatable[1..3,.]            // 取出参照决策单元的投入数据  
        Yref=datatable[4..5,.]            // 取出参照决策单元的产出数据  
        c = (1, J(1,N1,0))                // 目标函数的系数向量  
        lowerbd =., J(1,N1,0)             // 决策变量的下界  
        upperbd = J(1,N1+1,.)             // 决策变量上界  
        Aie = (J(3,1,0),Xref \ Y,-Yref)   // 不等式约束的系数矩阵  
        bie = X \ J(2,1,0)                // 不等式约束的上界  
        lp.setCoefficients(c)              // 设定lp的目标函数系数向量  
        lp.setInequality(Aie, bie)         // 设定lp的不等式约束  
        lp.setBounds(lowerbd, upperbd)     // 设定lp的决策变量上下界  
        theta=lp.optimize()                // 优化求解                 
        te=1/theta                        // 技术效率等于theta的倒数  
        te                                // 将te值列印出来         
end

*##12.4 步骤3：循环迭代

/////////

mata:   
       datatable=( 5,5,8,34,9     ///  
                 \ 6,9,6,67,6     ///  
                 \ 7,7,87,32,19   ///  
                 \ 23,14,45,11,31 ///  
                 \ 67,34,12,37,29 )  
	                    
	        lp = LinearProgram()             // 定义lp为线性规划类  
	        N1 = 5                           // 参照决策单元的个数  
	        N2 = 5                           // 评价决策单元的个数  
            theta=J(N2,1,.)                  // 定义theta用于储存计算结果  
        for(i=1;i<=N2;i++){                   
            X=datatable[1..3,i]              // 取出第i个决策单元的投入数据  
	        Y=datatable[4..5,i]              // 取出第i个决策单元的产出数据  
            Xref=datatable[1..3,.]           // 取出参照决策单元的投入数据  
	        Yref=datatable[4..5,.]           // 取出参照决策单元的产出数据  
            c = (1, J(1,N1,0))               // 目标函数的系数向量  
	        lowerbd =., J(1,N1,0)            //  决策变量的下界  
	        upperbd = J(1,N1+1,.)            //  决策变量上界  
            Aie = (J(3,1,0),Xref \ Y,-Yref)  // 不等式约束的系数矩阵  
            bie = X \ J(2,1,0)               // 不等式约束的上界  
            lp.setCoefficients(c)            // 设定lp的目标函数系数向量  
            lp.setInequality(Aie, bie)       // 设定lp的不等式约束  
	        lp.setBounds(lowerbd, upperbd)   // 设定lp的决策变量上下界  
            theta[i]=lp.optimize()           // 优化求解，将结果保存在theta
    }  
            te=1:/theta                      // 技术效率等于theta的倒数   
	        te                               // 将te值列印出来  
end  

/////////

mata:   
    real matrix function dea1(real matrix data, 
                              real matrix dataref, 
                              real scalar k       )  
   {  
        class LinearProgram scalar lp  
        lp = LinearProgram()             // 定义lp为线性规划类         
        N1 = cols(data)                  // 评价决策单元的个数  
        N2 = cols(dataref)               // 参照决策单元的个数  
        theta=J(N1,1,.)                     
        M=rows(data)                     // 投入产出变量个数和  
   for(i=1;i<=N1;i++){                   
        X=data[1..k,i]                   // 取出第i个决策单元的投入数据  
        Y=data[(k+1)..M,i]               // 取出第i个决策单元的产出数据  
        Xref=dataref[1..k,.]             // 取出参照决策单元的投入数据  
        Yref=dataref[(k+1)..M,.]         // 取出参照决策单元的产出数据  
        c = (1, J(1,N2,0))               // 目标函数的系数向量  
	    lowerbd =., J(1,N2,0)            // 决策变量的下界  
        upperbd = J(1,N2+1,.)            // 决策变量上界  
        Aie = (J(k,1,0),Xref \ Y,-Yref)  // 不等式约束的系数矩阵  
        bie = X \ J(M-k,1,0)             // 不等式约束的上界  
        lp.setCoefficients(c)            // 设定lp的目标函数系数向量  
        lp.setInequality(Aie, bie)       // 设定lp的不等式约束  
	    lp.setBounds(lowerbd, upperbd)   // 设定lp的决策变量上下界  
        theta[i]=lp.optimize()           // 优化求解，将结果保存在theta
	    }  
    te=1:/theta  
    return(te)  
	}  
end 

/////////

cap mata mata drop dea1()           //清除内存中可能存在的dea1()函数                                           
mata:
   real matrix function dea1(real matrix data, 
                             real matrix dataref, 
                              real scalar k       )  
    {  
        class LinearProgram scalar lp  
        lp = LinearProgram()             // 定义lp为线性规划类         
        N1 = cols(data)                  // 评价决策单元的个数  
        N2 = cols(dataref)               // 参照决策单元的个数  
	    theta=J(N1,1,.)                     
        M=rows(data)                     // 投入产出变量个数和  
	    for(i=1;i<=N1;i++){                   
        X=data[1..k,i]                   // 取出第i个决策单元的投入数据  
        Y=data[(k+1)..M,i]               // 取出第i个决策单元的产出数据  
        Xref=dataref[1..k,.]             // 取出参照决策单元的投入数据  
        Yref=dataref[(k+1)..M,.]         // 取出参照决策单元的产出数据  
        c = (1, J(1,N2,0))               // 目标函数的系数向量  
        lowerbd =., J(1,N2,0)            // 决策变量的下界  
        upperbd = J(1,N2+1,.)            // 决策变量上界  
        Aie = (J(k,1,0),Xref \ Y,-Yref)  // 不等式约束的系数矩阵  
	    bie = X \ J(M-k,1,0)             // 不等式约束的上界  
        lp.setCoefficients(c)            // 设定lp的目标函数系数向量  
        lp.setInequality(Aie, bie)       // 设定lp的不等式约束  
        lp.setBounds(lowerbd, upperbd)   // 设定lp的决策变量上下界  
        theta[i]=lp.optimize()           // 优化求解，将结果保存在theta
     }  
    te=1:/theta  
    return(te)  
}  
mata mosave dea1()                         // 将dea1()函数保存为本地文件
end 

/////////

mata:  
       datatable=(  5, 5, 8, 34, 9  ///  
                  \ 6, 9, 6, 67, 6  ///  
                  \ 7, 7, 87,32, 19 ///  
                  \ 23,14,45,11, 31 ///  
                  \ 67,34,12,37, 29    )  
    data=datatable  
    dataref=datatable  
    k=3  
    te=dea1(data,dataref,3)  
    te  
end

/////////

capture mata mata drop dea2()  
mata:   
    void function dea2(string scalar varname,
                       real   scalar k,
                       string scalar tename   )  
   {  
        data=st_data(.,varname)         // 将数据传递给Mata中的矩阵data  
        data=data'  
        dataref=data  
        class LinearProgram scalar lp  
	    lp = LinearProgram()            // 定义lp为线性规划类          
        N1=cols(data)                   // 评价决策单元的个数  
        N2=cols(dataref)                // 参照决策单元的个数  
        theta=J(N1,1,.)                     
        M=rows(data)                    // 投入产出变量个数和  
	   for(i=1;i<=N1;i++){                   
	    X=data[1..k,i]                  // 取出第i个决策单元的投入数据  
        Y=data[(k+1)..M,i]              // 取出第i个决策单元的产出数据  
        Xref=dataref[1..k,.]            // 取出参照决策单元的投入数据  
        Yref=dataref[(k+1)..M,.]        // 取出参照决策单元的产出数据  
        c = (1, J(1,N2,0))              // 目标函数的系数向量  
        lowerbd =., J(1,N2,0)           // 决策变量的下界  
        upperbd = J(1,N2+1,.)           // 决策变量上界  
        Aie = (J(k,1,0),Xref \ Y,-Yref) // 不等式约束的系数矩阵  
	    bie = X \ J(M-k,1,0)            // 不等式约束的上界  
	    lp.setCoefficients(c)           // 设定lp的目标函数系数向量  
        lp.setInequality(Aie, bie)      // 设定lp的不等式约束  
        lp.setBounds(lowerbd, upperbd)  // 设定lp的决策变量上下界  
        theta[i]=lp.optimize()          // 优化求解，将结果保存在theta
    }  
    st_view(te=.,.,tename)              // 将Stata变量映射到Mata中的te
    te[.,.]=1:/theta                    // 对te向量进行赋值，传递Stata 
	}  
    mata mosave dea2(),replace          // 将dea2()函数保存为本地文件  
end  

/////////

use Ex3.dta  
gen te=.  
mata:dea2("K L Y",2,"te") 

*##12.5 步骤4：将mata函数封装成Stata命令

capture program drop dea2       //清除内存中可能存在的dea2命令                                              
program define dea2  
    version 16  
    syntax, inputvars(varlist) outputvars(varlist) te(string)  
    confirm new var `te'  
    qui gen `te'=.  
    local nip: word count `inputvars'  
    mata:dea2("`inputvars' `outputvars'", `nip', "`te'")  
end  

/////////

use Ex3.dta  ,clear
dea2, inputvars(K L) outputvars(Y) te(te)                                     


/////////
program define dea3   
    version 16   
    syntax [if] [in], inputvars(varlist) outputvars(varlist) te(string)  
	marksample touse   
	confirm new var `te'   
	qui gen `te'=.   
    local nip: word count `inputvars'   
	mata:dea3("`inputvars' `outputvars'", "`touse'",`nip', "`te'")   
end   
cap mata mata drop dea3()   
mata:    
    void function dea3(string scalar varname,
                       string scalar flag,
	                       real   scalar k,
	                       string scalar tename)   
	    {   
	    data=st_data(.,varname,flag)    // 将Stata数据传递给Mata矩阵 
        data=data'   
        dataref=st_data(.,varname)   
        dataref=dataref'   
	    class LinearProgram scalar lp   
	    lp = LinearProgram()            // 定义lp为线性规划类           
        N1=cols(data)                   // 评价决策单元的个数   
        N2=cols(dataref)                // 参照决策单元的个数   
        theta=J(N1,1,.)                      
	    M=rows(data)                    // 投入产出变量个数和   
   for(i=1;i<=N1;i++){                    
	    X=data[1..k,i]                  // 取出第i个决策单元的投入数据   
        Y=data[(k+1)..M,i]              // 取出第i个决策单元的产出数据   
        Xref=dataref[1..k,.]            // 取出参照决策单元的投入数据   
        Yref=dataref[(k+1)..M,.]        // 取出参照决策单元的产出数据   
        c = (1, J(1,N2,0))              // 目标函数的系数向量   
        lowerbd =., J(1,N2,0)           // 决策变量的下界   
        upperbd = J(1,N2+1,.)           // 决策变量上界   
	    Aie = (J(k,1,0),Xref \ Y,-Yref) // 不等式约束的系数矩阵   
	    bie = X \ J(M-k,1,0)            // 不等式约束的上界   
        lp.setCoefficients(c)           // 设定lp的目标函数系数向量   
        lp.setInequality(Aie, bie)      // 设定lp的不等式约束   
        lp.setBounds(lowerbd, upperbd)  // 设定lp的决策变量上下界   
	    theta[i]=lp.optimize()          // 优化求解，将结果保存在theta
    }   
	 st_view(te=.,.,tename,flag)         // 将tename映射到Mata中的te   
	 te[.,.]=1:/theta                    // 对te进行赋值，传递回Stata 
}   
end  

/////////

use Ex3.dta  
dea3 if _n<3, inputvars(K L) outputvars(Y) te(te) 

/////////

	program define dea4  
    version 16  
	    syntax [if] [in], inputvars(varlist) outputvars(varlist) ///
	                      te(string) [reference(varname)]  
	    marksample touse  
	    confirm new var `te'  
	    if `"`reference'"'==""{  
	        tempvar reference  
	        qui gen byte `reference'=1        
	    }  
        else{  
        confirm numeric var `reference'  
        }  
    qui gen `te'=.  
    local nip: word count `inputvars'  
    mata:dea4("`inputvars' `outputvars'",  ///
	              "`touse'","`reference'",`nip', "`te'")  
	end  
cap mata mata drop dea4()  
mata:   
	    void function dea4(string scalar varname,
	                       string scalar flag,
                     string scalar rflag,
	                       real    scalar k,
                       string scalar tename)  
    {  
        data=st_data(.,varname,flag)    // 将Stata数据传递给Mata矩阵 
        data=data'  
        dataref=st_data(.,varname,rflag)  
	    dataref=dataref'  
        class LinearProgram scalar lp  
	    lp = LinearProgram()            // 定义lp为线性规划类          
	    N1 = cols(data)                 // 评价决策单元的个数  
	    N2 = cols(dataref)              // 参照决策单元的个数  
	    theta=J(N1,1,.)                     
        M=rows(data)                    // 投入产出变量个数和  
	   for(i=1;i<=N1;i++){                   
        X=data[1..k,i]                  // 取出第i个决策单元的投入数据  
	    Y=data[(k+1)..M,i]              // 取出第i个决策单元的产出数据  
        Xref=dataref[1..k,.]            // 取出参照决策单元的投入数据  
	    Yref=dataref[(k+1)..M,.]        // 取出参照决策单元的产出数据  
        c = (1, J(1,N2,0))              // 目标函数的系数向量  
	    lowerbd =., J(1,N2,0)           // 决策变量的下界  
	    upperbd = J(1,N2+1,.)           // 决策变量上界  
	    Aie = (J(k,1,0),Xref \ Y,-Yref) // 不等式约束的系数矩阵  
        bie = X \ J(M-k,1,0)            // 不等式约束的上界  
	    lp.setCoefficients(c)           // 设定lp的目标函数系数向量  
	    lp.setInequality(Aie, bie)      // 设定lp的不等式约束  
	    lp.setBounds(lowerbd, upperbd)  // 设定lp的决策变量上下界  
	    theta[i]=lp.optimize()          // 优化求解，将结果保存在theta
    }  
	st_view(te=.,.,tename,flag)         // 将tename映射到Mata中的te
    te[.,.]=1:/theta                    // 对te向量进行赋值，传递回Stata
}  
end

///////// 

use Ex3.dta  
gen rf=(_n<=10)  
dea4 if _n<3, inputvars(K L) outputvars(Y) te(te) reference(rf) 

/////////

local 0 K L = Y if _n<3, te(te) reference(rf)  
gettoken vars 0:0, p("=")  
di "`vars'"  
di "`0'"

/////////

gettoken vars 0:0, p("=")  
di "`vars'"  
di "`0'" 

/////////

program define dea5  
   version 16  
    gettoken inputvars 0:0, p("=")  
	    gettoken var 0:0, p("=")  
	    syntax varlist [if] [in], te(string) [reference(varname)]  
	    local outputvars `varlist'  
	    marksample touse  
	    confirm new var `te'  
    if `"`reference'"'==""{  
        tempvar reference  
	        qui gen byte `reference'=1        
    }  
	    else{  
	        confirm numeric var `reference'  
	    }  
	    qui gen `te'=.  
   local nip: word count `inputvars'  
	    mata:dea4("`inputvars' `outputvars'", ///
              "`touse'","`reference'",`nip', "`te'")  
end  
cap mata mata drop dea4()  
	mata:   
    void function dea4(string scalar varname,
	                   string scalar flag,   
	                   string scalar rflag,
	                   real   scalar k,
                       string scalar tename  )  
    {  
        data=st_data(.,varname,flag)       
        data=data'  
        dataref=st_data(.,varname,rflag)  
        dataref=dataref'  
        class LinearProgram scalar lp  
        lp = LinearProgram()                      
        N1 = cols(data)                      
        N2 = cols(dataref)                  
        theta=J(N1,1,.)                     
        M=rows(data)                      
   for(i=1;i<=N1;i++){                   
        X=data[1..k,i]                    
        Y=data[(k+1)..M,i]                 
        Xref=dataref[1..k,.]                
        Yref=dataref[(k+1)..M,.]           
	    c = (1, J(1,N2,0))                  
        lowerbd =., J(1,N2,0)              
        upperbd = J(1,N2+1,.)               
        Aie = (J(k,1,0),Xref \ Y,-Yref)     
	    bie = X \ J(M-k,1,0)               
        lp.setCoefficients(c)               
	    lp.setInequality(Aie, bie)          
        lp.setBounds(lowerbd, upperbd)      
        theta[i]=lp.optimize()               
    }  
    st_view(te=.,.,tename,flag)            
	te[.,.]=1:/theta                        
	}  
end      

/////////

use Ex3.dta ,clear 
gen rf=(_n<=10)  
dea5 K L=Y if _n<3, te(te) reference(rf) 

*#第十三章 DEA程序编写的实战演练
*##13.1 SBM模型的程序编写

	cap mata mata drop sbm()  
	mata:   
	    void function sbm(string scalar varname,
	                      string scalar flag,   
	                      string scalar rflag, 
	                      real   scalar p,   
	                      string scalar strname  )  
	    {  
	     data=st_data(.,varname,flag)       
	     data=data'  
	     dataref=st_data(.,varname,rflag)  
	     dataref=dataref'  
         class LinearProgram scalar lp  
         lp = LinearProgram()   
	     lp.setMaxOrMin("min")  
	        N1=cols(data)                       
	        N2=cols(dataref)                                      
	        M=rows(data)   
	        beta=J(N1,M+1,.)           
   for(i=1;i<=N1;i++){                   
	        X=data[1..p,i]                      
	        Y=data[(p+1)..M,i]                 
            Xref=dataref[1..p,.]                
	        Yref=dataref[(p+1)..M,.]            
	        c = (1, -1/p*(1:/X'),J(1,M-p,0),J(1,N2,0))                  
	        lowerbd = J(1,length(c),0)               
            upperbd = J(1,length(c),.)               
            Ae1 = (1,J(1,p,0),1/(M-p)*(1:/Y'),J(1,N2,0))     
	        Ae2 = (-X,I(p),J(p,M-p,0),Xref)     
	        Ae3 = (-Y,J(M-p,p,0),-I(M-p),Yref)               
	        lp.setCoefficients(c)                
	        lp.setEquality(Ae1\Ae2\Ae3, 1 \ J(M,1,0))           
	        lp.setBounds(lowerbd, upperbd)       
            theta=lp.optimize()  
	        if (lp.converged()==1){  
	            bslk=lp.parameters()  
	            beta[i,.]=theta,bslk[1,2..(M+1)]/bslk[1,1]          
	        }  
    }  
    st_view(teslack=.,.,strname,flag)              
	    teslack[.,.]=beta                         
}  
end  

/////////	
	
	program define sbm  
    version 16  
    gettoken inputvars 0:0, p("=")  
    gettoken var 0:0, p("=")  
	    syntax varlist [if] [in], te(string) slack(string) ///
	                              [reference(varname)]  
	    local outputvars `varlist'  
	    marksample touse  
	    confirm new var `te'  
	    if `"`reference'"'==""{  
	        tempvar reference  
	        qui gen byte `reference'=1        
    }  
	    else{  
	        confirm numeric `reference'  
    }  
	    qui gen `te'=.  
    foreach v in `inputvars' `outputvars'{  
        qui gen `slack'_`v'=.  
        label var `slack'_`v' "Slack in `v'"  
	        local slackvars `slackvars' `slack'_`v'  
    }  
	    local nip: word count `inputvars'  
	    mata:sbm("`inputvars' `outputvars'", ///
             "`touse'","`reference'",`nip', "`te' `slackvars'")  
	end  
  
	cap mata mata drop sbm()  
	mata:   
	    void function sbm(string scalar varname,string scalar flag,   
	                      string scalar rflag, real scalar p,   
	                      string scalar strname)  
	    {  
	     data=st_data(.,varname,flag)       
	     data=data'  
	     dataref=st_data(.,varname,rflag)  
		 dataref=dataref'  
		 class LinearProgram scalar lp  
		 lp = LinearProgram()   
	     lp.setMaxOrMin("min")  
	     N1=cols(data)                       
	     N2=cols(dataref)                                      
         M=rows(data)   
	     beta=J(N1,M+1,.)           
	   for(i=1;i<=N1;i++){                   
	        X=data[1..p,i]                      
	        Y=data[(p+1)..M,i]                 
	        Xref=dataref[1..p,.]                
			Yref=dataref[(p+1)..M,.]            
	        c = (1, -1/p*(1:/X'),J(1,M-p,0),J(1,N2,0))                  
			lowerbd = J(1,length(c),0)               
			upperbd = J(1,length(c),.)               
			Ae1 = (1,J(1,p,0),1/(M-p)*(1:/Y'),J(1,N2,0))     
			Ae2 = (-X,I(p),J(p,M-p,0),Xref)     
			Ae3 = (-Y,J(M-p,p,0),-I(M-p),Yref)               
	        lp.setCoefficients(c)                
	        lp.setEquality(Ae1\Ae2\Ae3, 1 \ J(M,1,0))           
	        lp.setBounds(lowerbd, upperbd)       
	        theta=lp.optimize()  
	        if (lp.converged()==1){  
	            bslk=lp.parameters()  
            beta[i,.]=theta,bslk[1,2..(M+1)]/bslk[1,1]          
        }  
    }  
    st_view(teslack=.,.,strname,flag)              
	    teslack[.,.]=beta                         
}  
	end 

/////////
	
use Ex3.dta ,clear 
sbm K L=Y, te(te) slack(S)

*##13.2 包含非合意产出的SBM模型的程序编写

	cap mata mata drop sbm2()  
	  mata:   
	      void function sbm2(string scalar varname,
                         string scalar flag,   
	                     string scalar rflag, 
                         real   scalar p,   
	                     real   scalar q,
                         string scalar strname )  
      {  
	       data=st_data(.,varname,flag)     
	       data=data'  
           dataref=st_data(.,varname,rflag)  
           dataref=dataref'  
           class LinearProgram scalar lp  
           lp = LinearProgram()   
           lp.setMaxOrMin("min")  
	       N1=cols(data)                       
	       N2=cols(dataref)                                      
           M=rows(data)   
	       beta=J(N1,M+1,.)           
	     for(i=1;i<=N1;i++){                 
          X=data[1..p,i]                      
          Y=data[(p+1)..(p+q),i]    
	      B=data[(p+q+1)..M,i]   
          Xref=dataref[1..p,.]                
          Yref=dataref[(p+1)..(p+q),.]    
          Bref=dataref[(p+q+1)..M,.]             
          c = (1, -1/p*(1:/X'),J(1,M-p,0),J(1,N2,0))                  
          lowerbd = J(1,length(c),0)               
          upperbd = J(1,length(c),.)               
          Ae1 = (1,J(1,p,0), 1/(M-p)*(1:/Y'), \\\
                 1/(M-p)*(1:/B'),J(1,N2,0)         )     
	      Ae2 = (-X,I(p),J(p,M-p,0),Xref)     
          Ae3 = (-Y,J(q,p,0),-I(q),J(q,M-p-q,0),Yref)  
          Ae4 = (-B,J(M-p-q,p+q,0),I(M-p-q),Bref)  
	          lp.setCoefficients(c)                
	          lp.setEquality(Ae1\Ae2\Ae3\Ae4, 1\J(M,1,0))           
	          lp.setBounds(lowerbd, upperbd)       
	          theta=lp.optimize()  
          if (lp.converged()==1){  
            bslk=lp.parameters()  
          beta[i,.]=theta,bslk[1,2..(M+1)]/bslk[1,1]          
          }  
     }  
	      st_view(teslack=.,.,strname,flag)            
    teslack[.,.]=beta                         
  }  
  end 

/////////
  
	program define sbm2  
    version 16  
    gettoken inputvars 0:0, parse("=")  
    gettoken var 0:0, parse("=")  
	    gettoken doutputvars 0:0, parse(":")  
    gettoken var 0:0, parse(":")  
	    syntax varlist [if] [in], te(string) slack(string) ///
                              [reference(varname)]  
	    local udoutputvars `varlist'  
    marksample touse  
    confirm new var `te'  
    if `"`reference'"'==""{  
	        tempvar reference  
        qui gen byte `reference'=1        
    }  
    else{  
	        confirm numeric `reference'  
	    }  
	    qui gen `te'=.  
    foreach v in `inputvars' `doutputvars' `udoutputvars'{  
	        qui gen `slack'_`v'=.  
	        label var `slack'_`v' "Slack in `v'"  
        local slackvars `slackvars' `slack'_`v'  
	    }  
	    local nx: word count `inputvars'  
	    local ny: word count `doutputvars'       
	    mata:sbm2("`inputvars' `doutputvars' `udoutputvars'", ///
	              "`touse'","`reference'",`nx',`ny', "`te' `slackvars'")  
	end  
	  
	cap mata mata drop sbm2()  
	mata:   
	    void function sbm2(string scalar varname,
                       string scalar flag,   
	                       string scalar rflag, 
	                       real   scalar p,   
	                       real   scalar q,
	                       string scalar strname  )  
    {  
	     data=st_data(.,varname,flag)       
     data=data'  
	     dataref=st_data(.,varname,rflag)  
	     dataref=dataref'  
	     class LinearProgram scalar lp  
	     lp = LinearProgram()   
	     lp.setMaxOrMin("min")  
	        N1=cols(data)                       
	        N2=cols(dataref)                                      
	        M=rows(data)   
	        beta=J(N1,M+1,.)           
	   for(i=1;i<=N1;i++){                   
	        X=data[1..p,i]                      
	        Y=data[(p+1)..(p+q),i]    
	        B=data[(p+q+1)..M,i]   
	        Xref=dataref[1..p,.]                
	        Yref=dataref[(p+1)..(p+q),.]    
	        Bref=dataref[(p+q+1)..M,.]             
        c = (1, -1/p*(1:/X'),J(1,M-p,0),J(1,N2,0))                  
	        lowerbd = J(1,length(c),0)               
	        upperbd = J(1,length(c),.)               
	        Ae1 = (1,J(1,p,0),1/(M-p)*(1:/Y'),1/(M-p)*(1:/B'),J(1,N2,0))     
	        Ae2 = (-X,I(p),J(p,M-p,0),Xref)     
	        Ae3 = (-Y,J(q,p,0),-I(q),J(q,M-p-q,0),Yref)  
	        Ae4 = (-B,J(M-p-q,p+q,0),I(M-p-q),Bref)  
	        lp.setCoefficients(c)                
        lp.setEquality(Ae1\Ae2\Ae3\Ae4, 1\J(M,1,0))           
        lp.setBounds(lowerbd, upperbd)       
        theta=lp.optimize()  
        if (lp.converged()==1){  
	            bslk=lp.parameters()  
	            beta[i,.]=theta,bslk[1,2..(M+1)]/bslk[1,1]          
	        }  
	    }  
	    st_view(teslack=.,.,strname,flag)              
	    teslack[.,.]=beta                         
	}  
	end 

/////////

use Ex4.dta ,clear 
sbm2 K L=Y:CO2, te(te) slack(S)

*##13.3 DDF模型的程序编写

cap mata mata drop ddf()  
   mata:   
	       void function ddf( string scalar varname,  
	                          string scalar gname,  
	                          string scalar flag,   
	                          string scalar rflag,   
	                          real   scalar p,  
	                          real   scalar q,  
	                          string scalar strname )  
       {  
	        data=st_data(.,varname,flag)        
	        data=data'  
			dataref=st_data(.,varname,rflag)  
			dataref=dataref'  
			gvec=st_data(.,gname,flag)  
			gvec=gvec'  
	        class LinearProgram scalar lp  
	        lp = LinearProgram()   
	        lp.setMaxOrMin("max")  
			N1=cols(data)                       
			N2=cols(dataref)                                      
			M=rows(data)   
	        beta=J(N1,1,.)    
	  
      for(i=1;i<=N1;i++){                    
	           X=data[1..p,i]                      
	           Y=data[(p+1)..(p+q),i]    
	           B=data[(p+q+1)..M,i]  
	           gX=gvec[1..p,i]  
               gY=gvec[(p+1)..(p+q),i]  
	           gB=gvec[(p+q+1)..M,i]  
	           Xref=dataref[1..p,.]                
	           Yref=dataref[(p+1)..(p+q),.]    
	           Bref=dataref[(p+q+1)..M,.]             
	           c = (1, J(1,N2,0))   
	           lp.setCoefficients(c)   
	           lowerbd = .,J(1,length(c)-1,0)               
		       upperbd = J(1,length(c),.)  
               lp.setBounds(lowerbd, upperbd)  
               Aie1 = (-gX,Xref)     
               Aie2 = (gY,-Yref)  
	           lp.setInequality(Aie1 \ Aie2, X \ -Y)                   
	           Ae = (-gB,Bref)               
               lp.setEquality(Ae, B)     
	           beta[i]=lp.optimize()  
	       }  
	       st_view(te=.,.,strname,flag)            
       te[.,.]=beta                       
	   }  
	   end   

/////////
   
program define ddf  
    version 16  
    gettoken inputvars 0:0, parse("=")  
	    gettoken var 0:0, parse("=")  
	    gettoken doutputvars 0:0, parse(":")  
	    gettoken var 0:0, parse(":")  
	    syntax varlist [if] [in], dv(string)  [gx(varlist) gy(varlist) ///
	                              gb(varlist) reference(varname)]  
	    local udoutputvars `varlist'  
	    marksample touse  
	    confirm new var `dv'  
    if `"`reference'"'==""{  
	        tempvar reference  
	        qui gen byte `reference'=1        
    }  
	    else{  
	        confirm numeric var `reference'  
	    }  
	    if `"`gx'"'==""{  
	        foreach v in `inputvars'{  
	            tempvar g`v'  
	            qui gen `g`v''=-`v'  
            local gx `gx' `g`v''  
	        }         
	    }  
  
	    if `"`gy'"'==""{  
        foreach v in `doutputvars'{  
            tempvar g`v'  
            qui gen `g`v''=`v'  
	            local gy `gy' `g`v''  
        }         
	    }         
	    if `"`gb'"'==""{  
        foreach v in `udoutputvars'{  
	            tempvar g`v'  
	            qui gen `g`v''=-`v'  
	            local gb `gb' `g`v''  
	        }         
	    }         
	    qui gen `dv'=.  
	    local nx: word count `inputvars'  
	    local ny: word count `doutputvars'        
	    mata:ddf("`inputvars' `doutputvars' `udoutputvars'", ///
	             "`gx' `gy' `gb'", "`touse'","`reference'", ///
	              `nx',`ny', "`dv'"                              )  
	end  
	  
  
	cap mata mata drop ddf()  
	   mata:   
       void function ddf( string scalar varname,  
                          string scalar gname,  
	                          string scalar flag,   
	                          string scalar rflag,   
                          real scalar   p,  
	                          real scalar   q,  
	                          string scalar strname)  
	       {  
	        data=st_data(.,varname,flag)        
	        data=data'  
	        dataref=st_data(.,varname,rflag)  
	        dataref=dataref'  
	        gvec=st_data(.,gname,flag)  
	        gvec=gvec'  
	        class LinearProgram scalar lp  
	        lp = LinearProgram()   
	        lp.setMaxOrMin("max")  
	        N1=cols(data)                       
	        N2=cols(dataref)                                      
        M=rows(data)   
        beta=J(N1,1,.)    

	      for(i=1;i<=N1;i++){                    
           X=data[1..p,i]                      
	           Y=data[(p+1)..(p+q),i]    
	           B=data[(p+q+1)..M,i]  
	           gX=gvec[1..p,i]  
	           gY=gvec[(p+1)..(p+q),i]  
	           gB=gvec[(p+q+1)..M,i]  
	           Xref=dataref[1..p,.]                
	           Yref=dataref[(p+1)..(p+q),.]    
	           Bref=dataref[(p+q+1)..M,.]             
	           c = (1, J(1,N2,0))   
	           lp.setCoefficients(c)   
			   lowerbd = .,J(1,length(c)-1,0)               
               upperbd = J(1,length(c),.)  
	           lp.setBounds(lowerbd, upperbd)  
	           Aie1 = (-gX,Xref)     
	           Aie2 = (gY,-Yref)  
	           lp.setInequality(Aie1 \ Aie2, X \ -Y)                   
           Ae = (-gB,Bref)               
           lp.setEquality(Ae, B)     
           beta[i]=lp.optimize()  
       }  
       st_view(te=.,.,strname,flag)            
	       te[.,.]=beta                       
	   }  
	   end    

/////////
	   
use Ex4.dta ,clear 
ddf K L=Y:CO2, dv(beta1)   
	gen gk=-1  
	gen gl=-1  
	gen gy=1  
	gen gb=-1  
ddf K L=Y:CO2, dv(beta2) gx(gk gl) gy(gy) gb(gb)  

*13.4 NDDF模型的程序编写

cap mata mata drop nddf()    
mata:     
 void function nddf( string scalar varname,    
	                     string scalar gname,   //方向向量  
	                     real   rowvector w,    // 权重向量  
	                     string scalar flag,     
	                     string scalar rflag,     
	                     real   scalar p,    
	                     real   scalar q,    
                     string scalar strname )    
	    {    
	        data=st_data(.,varname,flag)          
	        data=data'    
	        dataref=st_data(.,varname,rflag)    
	        dataref=dataref'    
	        gvec=st_data(.,gname,flag)    
	        gvec=gvec'    
	        class LinearProgram scalar lp    
	        lp = LinearProgram()     
	        lp.setMaxOrMin("max")    
	        N1=cols(data)                         
	        N2=cols(dataref)                                        
	        M=rows(data)     
			beta=J(N1,M+1,.)      
	        Xref=dataref[1..p,.]                  
	        Yref=dataref[(p+1)..(p+q),.]      
	        Bref=dataref[(p+q+1)..M,.]      
	      for(i=1;i<=N1;i++){                      
	           X=data[1..p,i]                        
	           Y=data[(p+1)..(p+q),i]      
	           B=data[(p+q+1)..M,i]    
	           gX=gvec[1..p,i]    
	           gY=gvec[(p+1)..(p+q),i]    
	           gB=gvec[(p+q+1)..M,i]               
	           c = (w, J(1,N2,0))     
	           lp.setCoefficients(c)         
	           lowerbd = J(1,length(c),0)  
	           upperbd = J(1,length(c),.)    
	           lp.setBounds(lowerbd, upperbd)    
	            Aie1=diag(-gX),J(p,rows(data)-p,0),Xref  
	            bie1=X  
	            Aie2=J(q,p,0),diag(gY),J(q,rows(data)-p-q,0),-Yref  
	            bie2=-Y  
	            lp.setInequality((Aie1 \ Aie2 ), (bie1 \ bie2)) 
	            Aec=J(rows(data)-p-q,p+q,0),diag(-gB),Bref  
	            bec=B   
	            lp.setEquality(Aec,bec)           
	            dv=lp.optimize()      
	           if(lp.converged()==1){  
	               paras = lp.parameters()  
	               beta[i,.]=dv,paras[1,1..M]  
            }  
	       }    
	       st_view(te=.,.,strname,flag)              
	       te[.,.]=beta                         
	   }    
	   end    

/////////	

	cap program drop nddf    
	program define nddf      
	    version 16      
	    gettoken inputvars 0:0, parse("=")      
	    gettoken var 0:0, parse("=")      
	    gettoken doutputvars 0:0, parse(":")      
	    gettoken var 0:0, parse(":")      
	    syntax varlist [if] [in], dv(string) beta(string)  ///    
	                              [wmat(name) gx(varlist)  ///  
	                               gy(varlist) gb(varlist) ///  
	                               reference(varname)]      
	    local udoutputvars `varlist'      
	    marksample touse      
	    confirm new var `dv'      
	    if `"`reference'"'==""{      
	        tempvar reference      
	        qui gen byte `reference'=1            
	    }      
	    else{      
	        confirm numeric var `reference'      
	    }      
	    if `"`gx'"'==""{      
	        foreach v in `inputvars'{      
	            tempvar g`v'      
	            qui gen `g`v''=-`v'      
	            local gx `gx' `g`v''      
	        }             
	    }      
	      
	    if `"`gy'"'==""{      
	        foreach v in `doutputvars'{      
	            tempvar g`v'      
	            qui gen `g`v''=`v'      
            local gy `gy' `g`v''      
	        }             
	    }             
	    if `"`gb'"'==""{      
	        foreach v in `udoutputvars'{      
	            tempvar g`v'      
	            qui gen `g`v''=-`v'      
            local gb `gb' `g`v''      
        }             
	    }      
    
    qui gen `dv'=.      
	    local nx: word count `inputvars'      
	    local ny: word count `doutputvars'     
    local nb : word count `udoutputvars'    
	   tempname wv    
	   if "`wmat'"==""{    
       mata: `wv'=J(1,`nx'+`nb'+`ny',1)    
   }     
	   else{    
	      mata: `wv'=st_matrix("`wmat'")     
  }    
    
 foreach v in `inputvars' `doutputvars' `udoutputvars'{      
        qui gen `beta'_`v'=.       
	        local slackvars `slackvars' `beta'_`v'      
	    }      
	 mata: nddf("`inputvars' `doutputvars' `udoutputvars'",  ///    
	            "`gx' `gy' `gb'", `wv', "`touse'",   ///    
            "`reference'",`nx',`ny', "`dv' `slackvars'" )      
	end  
	  
 	cap mata mata drop nddf()      
	mata:       
	 void function nddf( string scalar varname,      
	                     string scalar gname,   //方向向量    
					     real   rowvector w,    // 权重向量    
	                     string scalar flag,       
	                     string scalar rflag,       
	                     real   scalar p,      
	                     real   scalar q,      
	                     string scalar strname )      
    {      
	        data=st_data(.,varname,flag)            
			data=data'      
	        dataref=st_data(.,varname,rflag)      
	        dataref=dataref'      
	        gvec=st_data(.,gname,flag)      
			gvec=gvec'      
			class LinearProgram scalar lp      
	        lp = LinearProgram()       
	        lp.setMaxOrMin("max")      
	        N1=cols(data)                           
	        N2=cols(dataref)                                          
	        M=rows(data)       
	        beta=J(N1,M+1,.)        
	        Xref=dataref[1..p,.]                    
	        Yref=dataref[(p+1)..(p+q),.]        
	        Bref=dataref[(p+q+1)..M,.]        
	      for(i=1;i<=N1;i++){                        
		       X=data[1..p,i]                          
	           Y=data[(p+1)..(p+q),i]        
	           B=data[(p+q+1)..M,i]      
               gX=gvec[1..p,i]      
	           gY=gvec[(p+1)..(p+q),i]      
	           gB=gvec[(p+q+1)..M,i]                 
	           c = (w, J(1,N2,0))       
	           lp.setCoefficients(c)           
	           lowerbd = J(1,length(c),0)    
	           upperbd = J(1,length(c),.)      
	           lp.setBounds(lowerbd, upperbd)      
	            Aie1=diag(-gX),J(p,rows(data)-p,0),Xref    
	            bie1=X    
	            Aie2=J(q,p,0),diag(gY),J(q,rows(data)-p-q,0),-Yref    
	            bie2=-Y    
	            lp.setInequality((Aie1 \ Aie2 ), (bie1 \ bie2))   
	            Aec=J(rows(data)-p-q,p+q,0),diag(-gB),Bref    
	            bec=B     
	            lp.setEquality(Aec,bec)             
	            dv=lp.optimize()        
	           if(lp.converged()==1){    
	               paras = lp.parameters()    
	               beta[i,.]=dv,paras[1,1..M]    
	            }    
	       }      
	       st_view(te=.,.,strname,flag)                
	       te[.,.]=beta                           
	   }      
	   end      

/////////
	   
use Ex4.dta,clear   
matrix w = (1,1,1,1)     
nddf K L=Y:CO2, dv(dv) beta(b_) wmat(w)

*##13.5 曼奎斯特生产率指数的程序编写

	use Ex4.dta,clear  
	gen D11=.  // D^{t}(X_t,Y_t)
	gen flag=.  
	forvalues j=1/3{  
	    cap drop te  
	    replace flag=(t==`j')  
	    dea5 K L=Y if t==`j', te(te) reference(flag)  
	    qui replace D11=1/te if t==`j'  
	}  
	  
	gen D22=.  // D^{t+1}(X_{t+1},Y_{t+1})
	forvalues j=2/3{  
	    cap drop te  
	    replace flag=(t==`j')  
	    dea5 K L=Y if t==`j', te(te) reference(flag)  
	    qui replace D22=1/te if t==`j'  
	}  
	
	gen D12=.   // D^{t}(X_{t+1},Y_{t+1})
	forvalues j=1/3{  
	    cap drop te  
	    replace flag=(t==`j')  
	    dea5 K L=Y if t==`j'+1, te(te) reference(flag)  
	    qui replace D12=1/te if t==`j'+1  
	}  
	  
	gen D21=.  // D^{t+1}(X_{t},Y_{t})
	forvalues j=2/3{  
	    cap drop te  
	    replace flag=(t==`j')  
	    dea5 K L=Y if t==`j'-1, te(te) reference(flag)  
	    qui replace D21=1/te if t==`j'-1  
	}  
bys id (t): gen from=t[_n-1]  
	bys id (t): gen to=t  
	bys id (t): gen mpi=sqrt(D12/D11[_n-1]*D22/D21[_n-1])  
	sort id t  
	list id from to mpi if t>1, sep(0)

/////////	
	
	program define mpi  
	    version 16  
	    gettoken inputvars 0:0, p("=")  
	    gettoken var 0:0, p("=")  
	    syntax varlist [if] [in],  id(varname) t(varname) [sav(string)]  
	    local outputvars `varlist'  
	    marksample touse  
	    preserve  
	    qui keep `inputvars' `outputvars' `id' `t' `touse'  
	    tempvar D11 D21 D12 flag te t2  
	    qui gen `D11'=.  
	    qui gen `D12'=.  
	    qui gen `D21'=.  
	    qui gen `flag'=.  
	    qui egen `t2'=group(`t')  
	    qui su `t2',meanonly  
        local T=r(max)  
	      
	    * computing D_{t}(t) & D_{t+1}(t+1)  
	    forvalue j=1/`T'{  
	        qui replace `flag'=(`t2'==`j')  
	        cap drop `te'  
	        dea5 `inputvars'=`outputvars' if `t2'==`j' & `touse', ///
	                                      te(`te') reference(`flag')  
	        qui replace `D11'=1/`te' if `t2'==`j' & `touse'  
	    }  
	    * computing D_{t}(t+1)   
	    forvalue j=2/`T'{  
	        qui replace `flag'=(`t2'==`j'-1)  
	        cap drop `te'  
	        dea5 `inputvars'=`outputvars' if `t2'==`j' & `touse', ///
	                                      te(`te') reference(`flag')  
	        qui replace `D12'=1/`te' if `t2'==`j' & `touse'  
	    }     
	    * computing D_{t+1}(t)   
    forvalue j=2/`T'{  
	        qui replace `flag'=(`t2'==`j')  
	        cap drop `te'  
	        dea5 `inputvars'=`outputvars' if `t2'==`j'-1 & `touse', ///
	                                     te(`te') reference(`flag')  
	        qui replace `D21'=1/`te' if `t2'==`j'-1 & `touse'  
	    }     
	      
	    qui bys `id' (`t2'): gen from=`t'[_n-1]  
	    qui bys `id' (`t2'): gen to=`t'  
	    bys `id' (`t2'): gen mpi=sqrt(`D12'/`D11'[_n-1]*`D11'/`D21'[_n-1])  
    qui drop if `t2'==1 | `touse'==0  
	    sort `id' `t2'  
	    disp _n(2) "Total Factor Productivity Change"  
	    list `id' from to mpi, sep(0)   
	    di "Note: missing value indicates infeasible problem."      
	   if `"`sav'"'!=""{  
      label var mpi "Malmquist Productivity Index"  
	      qui keep `id' from to mpi  
	      save `sav'  
	   }  
      
    restore  
	      
	end 

/////////

	use Ex4.dta,clear  
	mpi K L=Y, id(id) t(t) sav(mpi1) 

/////////

    capture program drop mpi             //清除内存中可能存在的mpi命令
	program define mpi  
	    version 16  
		gettoken inputvars 0:0, p("=")  
	    gettoken var 0:0, p("=")  
	    syntax varlist [if] [in], id(varname) t(varname)  ///
	            [SEQuential WINdow(numlist integer max=1 >=1) sav(string)]  
	    local outputvars `varlist'  
	    marksample touse  
	    preserve  
	    qui keep `inputvars' `outputvars' `id' `t' `touse'    
	    tempvar D11 D21 D12 flag te t2  
	    qui gen `D11'=.  
		qui gen `D12'=.  
		qui gen `D21'=.  
		qui gen `flag'=.  
	    qui egen `t2'=group(`t')  
		qui su `t2',meanonly  
	    local T=r(max)  
	      
	    local relf=cond("`sequential"!="","<=","==")  
    if "`window'"!=""{  
        local band=(`window'-1)/2  
	    }  
	      
	    * computing D_{t}(t) & D_{t+1}(t+1)  
	    forvalue j=1/`T'{  
        if "`window'"==""{  
	           qui replace `flag'=(`t2' `relf' `j')  
        }  
        else{  
	           qui replace `flag'=(`t2'>=`j'-`window' & `t2'<=`j'+`window')  
	        }  
	          
        cap drop `te'  
	        dea5 `inputvars'=`outputvars' if `t2'==`j' & `touse', ///
                                        te(`te') reference(`flag')  
        qui replace `D11'=1/`te' if `t2'==`j' & `touse'  
    }  
    * computing D_{t}(t+1)   
    forvalue j=2/`T'{  
	          
       if "`window'"==""{  
            qui replace `flag'=(`t2' `relf' `j'-1)  
        }  
        else{  
            qui replace `flag'=(`t2'>=`j'-1-`window' &  ///
	                                `t2'<=`j'-1+`window'        )  
	        }         
          
        cap drop `te'  
        dea5 `inputvars'=`outputvars' if `t2'==`j' & `touse', ///
	                                         te(`te') reference(`flag')  
	        qui replace `D12'=1/`te' if `t2'==`j' & `touse'  
	    }     
	    * computing D_{t+1}(t)   
	    forvalue j=2/`T'{  
	  
        if "`window'"==""{  
            qui replace `flag'=(`t2' `relf' `j')  
	        }  
	        else{  
           qui replace `flag'=(`t2'>=`j'-`window' & `t2'<=`j'+`window')  
        }             
          
        cap drop `te'  
        dea5 `inputvars'=`outputvars' if `t2'==`j'-1 & `touse', ///
	                                          te(`te') reference(`flag')  
	        qui replace `D21'=1/`te' if `t2'==`j'-1 & `touse'  
	    }     
	      
	    qui bys `id' (`t2'): gen from=`t'[_n-1]  
	    qui bys `id' (`t2'): gen to=`t'  
    bys `id' (`t2'): gen mpi=sqrt(`D12'/`D11'[_n-1]*`D11'/`D21'[_n-1])  
	    qui drop if `t2'==1 & `touse'==0  
    sort `id' `t2'  
	    disp _n(2) "Total Factor Productivity Change"  
	    list `id' from to mpi, sep(0)   
	    di "Note: missing value indicates infeasible problem."      
	   if `"`sav'"'!=""{  
      label var mpi "Malmquist Productivity Index"  
	      qui keep `id' from to mpi  
	      save `sav'  
   }  
	      
    restore  
	      
	end

/////////

	use Ex4.dta,clear  
	gen D11=.  
	forvalues j=1/5{  
	    cap drop te  
	    dea5 K L=Y if t==`j', te(te)   
	    qui replace D11=1/te if t==`j'  
	}  
	bys id (t): gen from=t[_n-1]  
	bys id (t): gen to=t  
	bys id (t): gen mpi=D11/D11[_n-1]  
	sort id t  
	list id from to mpi if t>1, sep(0) 

/////////

	use Ex4.dta,clear  
	gen D11=.  
	gen flag=.  
	bys id (t): gen from=t[_n-1]  
	bys id (t): gen to=t  
	gen mpi=.  
	forvalues j=2/5{  
	    cap drop te  
	    qui replace flag=(t==`j' | t==`j'-1)  
	    dea5 K L=Y if t==`j'-1, te(te) reference(flag) 
		  
	    qui replace D11=1/te if t==`j'-1 
		cap drop te
	    dea5 K L=Y if t==`j', te(te) reference(flag)  
		  
    qui replace D11=1/te if t==`j'  
	cap drop te
	    qui bys id (t): replace mpi=D11/D11[_n-1]  
	}  
	  
	sort id t  
	list id from to mpi if t>1, sep(0) 
	
/////////

	cap program drop mpi  
	program define mpi  
	    version 16  
	    gettoken inputvars 0:0, p("=")  
	    gettoken var 0:0, p("=")  
	    syntax varlist [if] [in], id(varname) t(varname) [SEQuential ///  
	                              WINdow(numlist intege max=1 >=1)   ///  
	                              global BIennial sav(string)       ]  
	    local outputvars `varlist'  
	    marksample touse  
	    preserve  
	    qui keep `id' `t' `touse' `inputvars' `outputvars'  
	    tempvar D11 D21 D12 flag te t2  
	    qui gen `D11'=.  
	    qui gen `D12'=.  
	    qui gen `D21'=.  
	    qui gen `flag'=.  
	    qui egen `t2'=group(`t')  
	    qui su `t2',meanonly  
	    local T=r(max)  
	    qui bys `id' (`t2'): gen from=`t'[_n-1]  
	    qui bys `id' (`t2'): gen to=`t'  
	  
	    if "`global'"!=""{  
	        forvalues j=1/`T'{  
	            cap drop `te'  
            qui dea5 `inputvars'=`outputvars' if `t2'==`j' & `touse', ///
	                                                  te(`te')   
	            qui replace `D11'=1/`te' if `t2'==`j' & `touse'  
	        }  
	        qui bys `id' (`t2'): gen mpi=`D11'/`D11'[_n-1]    
	    }  
	    else if "`biennial'"!=""{  
	        qui gen mpi=.  
	        forvalues j=2/`T'{  
	          cap drop `te'  
	          qui replace `flag'=(`t2'==`j' | `t2'==`j'-1)  
	          qui dea5 `inputvars'=`outputvars' if `t2'==`j'-1 & `touse', ///
	                                           te(`te') reference(`flag')  
	          qui replace `D11'=1/`te' if `t2'==`j'-1 & `touse'
			  cap drop `te'  
	          qui dea5 `inputvars'=`outputvars' if `t2'==`j' & `touse', ///
	                                          te(`te') reference(`flag')  
	          qui replace `D11'=1/`te' if `t2'==`j' & `touse'  
			  cap drop `te'   
	          qui bys `id' (`t2'): replace mpi=`D11'/`D11'[_n-1] ///
	                                          if `t2'==`j'       
	        }  
	          
	    }  
	    else{  
	        local relf=cond("`sequential'"!="","<=","==")  
	  
	        * computing D_{t}(t) & D_{t+1}(t+1)  
	        forvalue j=1/`T'{  
	            if "`window'"==""{  
	                qui replace `flag'=(`t2' `relf' `j')  
	            }  
	            else{  
	                qui replace `flag'=(`t2'>=`j'-`window' &  /// 
	                                          `t2'<=`j'+`window')  
	            }  
	  
	        cap drop `te'  
            dea5 `inputvars'=`outputvars' if `t2'==`j' & `touse',  ///
	                                          te(`te') reference(`flag')  
	            qui replace `D11'=1/`te' if `t2'==`j' & `touse'  
	        }  
	        * computing D_{t}(t+1)   
        forvalue j=2/`T'{  
	  
            if "`window'"==""{  
	                qui replace `flag'=(`t2' `relf' `j'-1)  
            }  
	            else{  
	                qui replace `flag'=(`t2'>=`j'-1-`window' &  ///
	                                            `t2'<=`j'-1+`window')  
	            }         
	
	            cap drop `te'  
	            dea5 `inputvars'=`outputvars' if `t2'==`j' & `touse',  ///
	                                            te(`te') reference(`flag')  
	            qui replace `D12'=1/`te' if `t2'==`j' & `touse'  
	        }     
	        * computing D_{t+1}(t)   
	        forvalue j=2/`T'{  
	  
	            if "`window'"==""{  
                qui replace `flag'=(`t2' `relf' `j')  
	            }  
	            else{  
	                qui replace `flag'=(`t2'>=`j'-`window' &  ///
	                                            `t2'<=`j'+`window')  
	            }             
	  
	            cap drop `te'  
	            dea5 `inputvars'=`outputvars' if `t2'==`j'-1 & `touse',  ///
	                                            te(`te') reference(`flag')  
	            qui replace `D21'=1/`te' if `t2'==`j'-1 & `touse'  
	        }     
	    bys `id' (`t2'): gen mpi=sqrt(`D12'/`D11'[_n-1]*`D11'/`D21'[_n-1])      
	    }  
	    qui drop if `t2'==1 & `touse'==0  
	    sort `id' `t2'  
	    disp _n(2) "Total Factor Productivity Change"  
	    list `id' from to mpi, sep(0)   
	    di "Note: missing value indicates infeasible problem."      
	   if `"`sav'"'!=""{  
	      label var mpi "Malmquist Productivity Index"  
	      qui keep `id' from to mpi  
	      save `sav'  
	   }  
	      
	    restore  
	      
	end  

/////////

use Ex4.dta,clear  
  
	* 计算时序技术下的mpi  
	mpi K L=Y, id(id) t(t) seq  
	  
	* 计算window(3)技术下的mpi  
	mpi K L=Y, id(id) t(t) window(3)  
	  
	* 计算全局技术下的mpi  
	mpi K L=Y, id(id) t(t) global  
	  
	* 计算两期技术下的mpi  
	mpi K L=Y, id(id) t(t) bi 

*13.6 曼奎斯特-卢恩博格生产率指数的编写

	use Ex4.dta,clear   
	  
	gen D11=.  
	gen flag=.  
	  
	forvalues j=1/4{   
	    cap drop dv  
	    replace flag=(t==`j')  
	    ddf K L=Y:CO2 if t==`j', dv(dv) reference(flag)   
	    qui replace D11=dv if t==`j'  
	}  
	  
	gen D22=.  
	forvalues j=2/5{    
	    cap drop dv  
	    replace flag=(t==`j')  
	    ddf K L=Y:CO2 if t==`j', dv(dv) reference(flag)   
	    qui replace D22=dv if t==`j'  
	}  
	  
	gen D12=.  
	forvalues j=1/4{    
	    cap drop dv  
	    replace flag=(t==`j')  
	    ddf K L=Y:CO2 if t==`j'+1, dv(dv) reference(flag)   
	    qui replace D12=dv if t==`j'+1  
}  
	  
	gen D21=.  
	forvalues j=2/5{    
	    cap drop dv  
	    replace flag=(t==`j')  
    ddf K L=Y:CO2 if t==`j'-1, dv(dv) reference(flag)   
	    qui replace D21=dv if t==`j'-1  
	}  
	bys id (t): gen from=t[_n-1]  
	bys id (t): gen to=t  
	bys id (t): gen mlpi=sqrt((1+D11[_n-1])/(1+D12)*(1+D21[_n-1])/(1+D22))  
	sort id t  
	list id from to mlpi if t>1, sep(0) 

/////////
	
	cap program drop mlpi  
	program define mlpi  
	    version 16  
	    gettoken inputvars 0:0, parse("=")  
	    gettoken var 0:0, parse("=")  
	    gettoken doutputvars 0:0, parse(":")  
	    gettoken var 0:0, parse(":")  
	    syntax varlist [if] [in], id(varname) t(varname) [SEQuential ///  
	                              WINdow(numlist integer max=1 >=1)   ///  
	                              global BIennial sav(string)       ]  
	    local udoutputvars `varlist'  
	    marksample touse  
	    preserve  
	    keep `id' `t' `inputvars' `doutputvars' ///
	         `udoutputvars' `gx' `gy' `gb' `touse'  
	    tempvar D11 D21 D12 flag te t2  
	    qui gen `D11'=.  
	    qui gen `D12'=.  
	    qui gen `D21'=.  
	    qui gen `flag'=.  
	    qui egen `t2'=group(`t')  
	    qui su `t2',meanonly  
	    local T=r(max)  
	    qui bys `id' (`t2'): gen from=`t'[_n-1]  
	    qui bys `id' (`t2'): gen to=`t'  
	  
	    if "`global'"!=""{  
	        forvalues j=1/`T'{  
	            cap drop `te'  
	            qui ddf `inputvars'=`doutputvars' : `udoutputvars' ///  
	                         if `t2'==`j' & `touse', dv(`te')         
	            qui replace `D11'=1/`te' if `t2'==`j' & `touse'  
	        }  
	        qui bys `id' (`t2'): gen mlpi=(`D11'[_n-1]+1)/(`D11'+1)  
	    }  
	    else if "`biennial'"!=""{  
	        qui gen mlpi=.  
	        forvalues j=2/`T'{  
	            cap drop `te'  
            qui replace `flag'=(`t2'==`j' | `t2'==`j'-1)  
	            qui ddf `inputvars'=`doutputvars' : `udoutputvars' ///  
	                if `t2'==`j'-1 & `touse', dv(`te') reference(`flag')  
	            qui replace `D11'=`te' if `t2'==`j'-1 & `touse'  
	            cap drop `te'  
				qui ddf `inputvars'=`doutputvars' : `udoutputvars' ///  
	                 if `t2'==`j' & `touse', dv(`te') reference(`flag')  
	            qui replace `D11'=`te' if `t2'==`j' & `touse'  
	            cap drop `te'  
				qui bys `id' (`t2'): replace  mlpi= ///
	                          (`D11'[_n-1]+1)/(`D11'+1) if `t2'==`j'         
	        }  
	    }  
	    else{  
	        local relf=cond("`sequential'"!="","<=","==")  
	        * computing D_{t}(t) & D_{t+1}(t+1)  
	        forvalue j=1/`T'{  
	          if "`window'"==""{  
	           qui replace `flag'=(`t2' `relf' `j')  
	          }  
	          else{  
	           qui replace `flag'=(`t2'>=`j'-`window' & `t2'<=`j'+`window')  
	          }  
	  
	          cap drop `te'  
          qui ddf `inputvars'=`doutputvars' : `udoutputvars' ///  
	                    if `t2'==`j' & `touse', dv(`te')  reference(`flag')  
	          qui replace `D11'=`te' if `t2'==`j' & `touse'  
	        }  
	        * computing D_{t}(t+1)   
	        forvalue j=2/`T'{  
	  
           if "`window'"==""{  
	            qui replace `flag'=(`t2' `relf' `j'-1)  
	            }  
	           else{  
	            qui replace `flag'=(`t2'>=`j'-1-`window' & ///
	                                `t2'<=`j'-1+`window'       )  
	            }         
	  
	            cap drop `te'  
            qui ddf `inputvars'=`doutputvars' : `udoutputvars' ///  
                if `t2'==`j' & `touse', dv(`te')  reference(`flag')  
            qui replace `D12'=`te' if `t2'==`j' & `touse'  
        }     
	        * computing D_{t+1}(t)   
        forvalue j=2/`T'{  
  
            if "`window'"==""{  
	              qui replace `flag'=(`t2' `relf' `j')  
	            }  
	            else{  
              qui replace `flag'=(`t2'>=`j'-`window' & ///
	                                  `t2'<=`j'+`window'       )  
	            }             
            cap drop `te'  
	            qui ddf `inputvars'=`doutputvars' : `udoutputvars' ///  
	                if `t2'==`j'-1 & `touse', dv(`te') reference(`flag')  
	            qui replace `D21'=`te' if `t2'==`j'-1 & `touse'  
	        }     
	     
	        qui bys `id' (`t2'): gen mlpi= ///
	         sqrt((1+`D11'[_n-1])/(1+`D12')*(1+`D21'[_n-1])/(1+`D11'))      
	    }  
	      
	  
	    qui drop if `t2'==1 | `touse'==0  
	    sort `id' `t2'  
	    disp _n(2) " Total Factor Productivity Change"  
	    list `id' from to mlpi, sep(0)   
	    di "Note: missing value indicates infeasible problem."      
	   if `"`sav'"'!=""{  
	      label var mlpi "Malmquist-Luenberger Productivity Index"  
	      qui keep `id' from to mpi  
	      save `sav'  
	   }  
	      
	    restore  
	      
	end 

*#第十四章 基于DEA对偶问题的影子价格计算
*##14.3 基于SBM对偶问题的影子价格计算程序
	
	use Ex4.dta,clear  
	putmata data = (K L Y CO2)  ,replace
	mata:  
	        data = data'  
	        P=2  
	        Q=1  
	        R=1  
	        s=Q+R  
	        Xref = data[1..P,.]  
	        Yref = data[(P+1)..(P+Q),.]  
	        Bref = data[(P+Q+1)..(P+Q+R),.]  
	        X    = data[1..P,1]  
	        Y    = data[(P+1)..(P+Q),1]   
	        B    = data[(P+Q+1)..(P+Q+R),1]          
	        c=(Y',-X',-B')  
	        lowerbd=J(1,length(c),0)   
	        upperbd=J(1,length(c),.)  
	        lp = LinearProgram()  
	        A1 = (Yref',-Xref',-Bref')  
	        b1 = J(cols(Xref),1,0)  
	        A2 = (J(P,Q,0),-P*I(P),J(P,R,0))  
	        b2 =-1:/X  
	        A3 = (-s*diag(Y)+J(Q,1,Y'),J(Q,1,-X'),J(Q,1,-B'))  
	        b3 =-J(Q,1,1)  
	        A4 = (J(R,1,Y'),J(R,1,-X'),-s*diag(B)+J(R,1,-B'))  
	        b4 =-J(R,1,1)                     
	        Aie=A1 \ A2 \ A3 \ A4  
	        bie=b1 \ b2 \ b3 \ b4  
	        lp.setCoefficients(c)  
	        lp.setInequality(Aie, bie)  
	        lp.setBounds(lowerbd, upperbd)  
	        z=lp.optimize()  //求解线性规划问题  
	        lp.parameters()  //显示最优解  
	end  

/////////

	use Ex4.dta,clear  
	putmata data = (K L Y CO2)  ,replace
	 mata:  
	        data = data'  
	        P=2  
	        Q=1  
	        R=1  
	        s=Q+R  
	        Xref = data[1..P,.]  
	        Yref = data[(P+1)..(P+Q),.]  
	        Bref = data[(P+Q+1)..(P+Q+R),.]  
	        beta = J(cols(data), P+Q+R,.) // 初始化矩阵beta，用以储存最优解
	        lp   = LinearProgram()  			
	        for(i=1;i<=cols(data);i++){  
	            X    = data[1..P,i]  
	            Y    = data[(P+1)..(P+Q),i]  
	            B    = data[(P+Q+1)..(P+Q+R),i]        
	            c=(Y',-X',-B')  
	            lowerbd=J(1,length(c),0)  
	            upperbd=J(1,length(c),.)  
	            A1 = (Yref',-Xref',-Bref')  
	            b1 = J(cols(Xref),1,0)  
	            A2 = (J(P,Q,0),-P*I(P),J(P,R,0))  
				b2 =-1:/X  
	            A3 = (-s*diag(Y)+J(Q,1,Y'),J(Q,1,-X'),J(Q,1,-B'))  
	            b3 =-J(Q,1,1)  
	            A4 = (J(R,1,Y'),J(R,1,-X'),-s*diag(B)+J(R,1,-B'))  
				b4 =-J(R,1,1)                     
	            Aie=A1 \ A2 \ A3 \ A4  
	            bie=b1 \ b2 \ b3 \ b4  
	            lp.setCoefficients(c)                                          
	            lp.setInequality(Aie, bie)  
	            lp.setBounds(lowerbd, upperbd)  
            z=lp.optimize()  
            if(lp.converged()==1){ //判断线性规划问题是否求得最优解  
	                beta[i,.]=lp.parameters()  
	              }  
	          }  
	            beta //显示所有求解结果  
	  end 
	  
/////////
	  
mata:  
void function sbmdual(string scalar vars,  
                      real   scalar P,  
                      real   scalar Q,  
                      string scalar gen)  
{  
    data= st_data(.,vars)  
    data=data'  
    dataref= st_data(.,vars)  
    dataref=dataref'      
    R=rows(data)-P-Q  
    s=Q+R      
        Xref = dataref[1..P,.]  
        Yref = dataref[(P+1)..(P+Q),.]  
        Bref = dataref[(P+Q+1)..(P+Q+R),.]  
        class LinearProgram scalar lp  
        beta = J(cols(data), P+Q+R,.)  
         lp   = LinearProgram() 
        for(i=1;i<=cols(data);i++){  
            X    = data[1..P,i]  
            Y    = data[(P+1)..(P+Q),i]  
            B    = data[(P+Q+1)..(P+Q+R),i]   
            c=(Y',-X',-B')  
            lowerbd=J(1,length(c),0)  
            upperbd=J(1,length(c),.)   
            A1 = (Yref',-Xref',-Bref')  
            b1 = J(cols(Xref),1,0)  
            A2 = (J(P,Q,0),-P*I(P),J(P,R,0))  
            b2 =-1:/X  
            A3 = (-s*diag(Y)+J(Q,1,Y'),J(Q,1,-X'),J(Q,1,-B'))  
            b3 =-J(Q,1,1)  
            A4 = (J(R,1,Y'),J(R,1,-X'),-s*diag(B)+J(R,1,-B'))  
            b4 =-J(R,1,1)                     
            Aie=A1 \ A2 \ A3 \ A4  
            bie=b1 \ b2 \ b3 \ b4  
            lp.setCoefficients(c)  
            lp.setInequality(Aie, bie)  
            lp.setBounds(lowerbd, upperbd)  
            theta=lp.optimize()  
            if(lp.converged()==1){  
                beta[i,.]=lp.parameters()  
              }  
            }  
            st_view(res=.,.,gen)  
            res[.,.]=beta  
}  
end 


	 
/////////

	gen pk=.  
	gen pl=.  
	gen py=.  
	gen pco2=.  
	mata: sbmdual("K L Y CO2",2,1,"py pk pl pco2")                             
	br  

/////////
cap mata mata drop sbndual()	
mata:  
void function sbmdual(string scalar vars,  
	                  string scalar touse,  
	                  string scalar rflag,  
                      real   scalar P,  
                      real   scalar Q,  
                      string scalar gen)  
{  
    data= st_data(.,vars,touse)  
    data=data'  
    dataref= st_data(.,vars,rflag)  
    dataref=dataref'      
    R=rows(data)-P-Q  
    s=Q+R      
        Xref = dataref[1..P,.]  
        Yref = dataref[(P+1)..(P+Q),.]  
        Bref = dataref[(P+Q+1)..(P+Q+R),.]  
        class LinearProgram scalar lp  
        beta = J(cols(data), P+Q+R,.)  
         lp   = LinearProgram() 
        for(i=1;i<=cols(data);i++){  
            X    = data[1..P,i]  
            Y    = data[(P+1)..(P+Q),i]  
            B    = data[(P+Q+1)..(P+Q+R),i]   
            c=(Y',-X',-B')  
            lowerbd=J(1,length(c),0)  
            upperbd=J(1,length(c),.)   
            A1 = (Yref',-Xref',-Bref')  
            b1 = J(cols(Xref),1,0)  
            A2 = (J(P,Q,0),-P*I(P),J(P,R,0))  
            b2 =-1:/X  
            A3 = (-s*diag(Y)+J(Q,1,Y'),J(Q,1,-X'),J(Q,1,-B'))  
            b3 =-J(Q,1,1)  
            A4 = (J(R,1,Y'),J(R,1,-X'),-s*diag(B)+J(R,1,-B'))  
            b4 =-J(R,1,1)                     
            Aie=A1 \ A2 \ A3 \ A4  
            bie=b1 \ b2 \ b3 \ b4  
            lp.setCoefficients(c)  
            lp.setInequality(Aie, bie)  
            lp.setBounds(lowerbd, upperbd)  
            theta=lp.optimize()  
            if(lp.converged()==1){  
                beta[i,.]=lp.parameters()  
              }  
            }  
            st_view(res=.,.,gen,touse)  
            res[.,.]=beta  
}  
end 


/////////

capture program drop sbmdual  
	program define sbmdual  
	    version 16  
	       gettoken invars 0:0, parse("=")  
	       gettoken var 0:0, parse("=")  
	       gettoken gopvars 0:0, parse(":")  
	       gettoken var 0:0, parse(":")  
	       syntax varlist [if] [in],gen(string) [rflag(varname)]     
	       local bopvars `varlist'  
	        marksample touse   
	        markout `touse' `invars' `gopvars' `bopvars'  
	        if `"`rflag'"'==""{  
	            tempvar rflag  
	            qui gen byte `rflag' = 1  
	        }  
	        tempvar touse2  
	        mark `touse2' if `rflag'  
	        markout `touse2' `invars' `gopvars' `bopvars'     
	        local data `invars' `gopvars' `bopvars'  
	        foreach v in `gopvars' `invars' `bopvars'{  
	            qui gen double `gen'`v' = .  
	            local gens `gens' `gen'`v'  
	        }  
	        local nx: word count `invars'  
	        local ny: word count `gopvars'  
	        mata: sbmdual("`data'","`touse'", "`touse2'",`nx',`ny',"`gens'")  
	  
	end     
	  
	cap mata mata drop sbmdual()  
	mata:  
	void function sbmdual(string scalar vars,  
	                      string scalar touse,  
	                      string scalar rflag,  
	                      real   scalar P,  
	                      real   scalar Q,  
	                      string scalar gen)  
	{  
	    data= st_data(.,vars,touse)  
	    data=data'  
	    dataref= st_data(.,vars,rflag)  
	    dataref=dataref'      
	    R=rows(data)-P-Q  
	    s=P+Q      
	        Xref = dataref[1..P,.]  
        Yref = dataref[(P+1)..(P+Q),.]  
        Bref = dataref[(P+Q+1)..(P+Q+R),.]  
	        class LinearProgram scalar lp  
	        beta = J(cols(data), P+Q+R,.)  
	        for(i=1;i<=cols(data);i++){  
	            X    = data[1..P,i]  
	            Y    = data[(P+1)..(P+Q),i]  
	            B    = data[(P+Q+1)..(P+Q+R),i]   
	            c=(Y',-X',-B')  
	            lowerbd=J(1,length(c),0)  
	            upperbd=J(1,length(c),.)  
	            lp = LinearProgram()  
	            A1 = (Yref',-Xref',-Bref')  
	            b1 = J(cols(Xref),1,0)  
            A2 = (J(P,Q,0),-P*I(P),J(P,R,0))  
	            b2 =-1:/X  
	            A3 = (-s*diag(Y)+J(Q,1,Y'),J(Q,1,-X'),J(Q,1,-B'))  
            b3 =-J(Q,1,1)  
	            A4 = (J(R,1,Y'),J(R,1,-X'),-s*diag(B)+J(R,1,-B'))  
	            b4 =-J(R,1,1)                     
	            Aie=A1 \ A2 \ A3 \ A4  
            bie=b1 \ b2 \ b3 \ b4  
	            lp.setCoefficients(c)   
	            lp.setInequality(Aie, bie)  
	            lp.setBounds(lowerbd, upperbd)  
            theta=lp.optimize()  
	            if(lp.converged()==1){  
	                beta[i,.]=lp.parameters()  
	              }  
	            }  
	            st_view(res=.,.,gen,touse)  
	            res[.,.]=beta  
	}  
	  
	end  

/////////

	use Ex4.dta,clear  
	sbmdual K L = Y:CO2, gen(P_)                                               //报错
    /*
	    sbmdual():  3200  conformability error
                 <istmt>:     -  function returned error
	*/
	
*##14.5 基于NDDF对偶问题的影子价格计算程序
	
	use Ex5.dta,clear    
	putmata data = (labor capital energy gdp co2),replace  
	mata:    
	        w=J(cols(data),1,1)  //设定权重都为1 
	        P=3  
	        Q=1  
	        R=1  
	        data =data'  
	        Xref = data[1..P,.]  
	        Yref = data[(P+1)..(P+Q),.]  
	        Bref = data[(P+Q+1)..(P+Q+R),.]  
			X=data[1..P,1]               //决策单元1的投入数据
	        Y=data[(P+1)..(P+Q),1]      //决策单元1的合意产出数据
			B=data[(P+Q+1)..(P+Q+R),1]  //决策单元1的非合意产出数据        
	        gx=-X  
	        gy=Y  
			gb=-B  
	        wx=w[1::P]  
	        wy=w[(P+1)::(P+Q)]  
	        wb=w[(P+Q+1)::(P+Q+R)]  
	          
	        c=(X',-Y',B')  
	        lowerbd=J(1,length(c),0)  
	        upperbd=J(1,length(c),.)  
	        lp = LinearProgram()  
	        lp.setMaxOrMin("min")  
	         
	        A1 = (-Xref',Yref',-Bref')  
	        b1 = J(cols(Xref),1,0)  
	        A2 = (diag(gx),J(P,Q,0),J(P,R,0))  
	          
	        b2 =-wx  
	        A3 = (J(Q,P,0),-diag(gy),J(Q,R,0))  
			b3 =-wy  
	          
	        A4 = (J(R,P,0),J(R,Q,0),diag(gb))  
			b4 =-wb                     
	        Aie=A1 \ A2 \A3 \A4  
	        bie=b1 \ b2 \b3 \ b4  
			lp.setCoefficients(c)  
	        lp.setInequality(Aie, bie)  
	        lp.setBounds(lowerbd, upperbd)  
	        lp.optimize()  
	        lp.parameters()  
end 

/////////

    * 对所有决策单位的对偶问题进行求解
	use Ex5.dta,clear      
	putmata data = (labor capital energy gdp co2),replace    
	mata:      
	        w=J(cols(data),1,1)  //设定权重都为1   
	        P=3    
	        Q=1    
	        R=1    
	        data =data'    
	        Xref = data[1..P,.]    
	        Yref = data[(P+1)..(P+Q),.]    
	        Bref = data[(P+Q+1)..(P+Q+R),.]    
	        beta = J(cols(data),rows(data),.)  
            lp =  LinearProgram()    
	        lp.setMaxOrMin("min") 		
	    for(i=1;i<=cols(data);i++){
	        X=data[1..P,i]               
	        Y=data[(P+1)..(P+Q),i]       
	        B=data[(P+Q+1)..(P+Q+R),i]           
	        gx=-X    
	        gy=Y    
	        gb=-B    
	        wx=w[1::P]    
	        wy=w[(P+1)::(P+Q)]    
	        wb=w[(P+Q+1)::(P+Q+R)]    	            
	        c=(X',-Y',B')    
	        lowerbd=J(1,length(c),0)    
	        upperbd=J(1,length(c),.)                                            	   
	        A1 = (-Xref',Yref',-Bref')    
	        b1 = J(cols(Xref),1,0)    
			A2 = (diag(gx),J(P,Q,0),J(P,R,0))    	            
	        b2 =-wx    
	        A3 = (J(Q,P,0),-diag(gy),J(Q,R,0))    
			b3 =-wy               
	        A4 = (J(R,P,0),J(R,Q,0),diag(gb))    
	        b4 =-wb                       
			Aie=A1 \ A2 \A3 \A4    
	        bie=b1 \ b2 \b3 \ b4    
			lp.setCoefficients(c)    
	        lp.setInequality(Aie, bie)    
			lp.setBounds(lowerbd, upperbd)    
	        dv=lp.optimize()    
	        if(lp.converged()==1){  
	          beta[i,.]= lp.parameters()  
	       }  
	    }    
	end 

/////////

cap mata mata drop nddfdual()  
   mata:   
       void function nddfdual( string scalar varname, 
                          string scalar flag,   
                          string scalar rflag,   
                          real   scalar P,  
                          real   scalar Q,  
                          real  colvector w,
                          string scalar gname,  //方向向量
                          string scalar strname )  
       {  
        data=st_data(.,varname,flag)        
        data=data'  
        dataref=st_data(.,varname,rflag)  
        dataref=dataref'  
        M=rows(data)
        R=M-P-Q
        Xref = dataref[1..P,.]
        Yref = dataref[(P+1)..(P+Q),.]
        Bref = dataref[(P+Q+1)..(P+Q+R),.]            
        gvec=st_data(.,gname,flag)  
        gvec=gvec'  
        class LinearProgram scalar lp  
        lp = LinearProgram()   
        lp.setMaxOrMin("min")  
		wx   = w[1..P,1]
		wy   = w[(P+1)..(P+Q),1]
		wb   = w[(P+Q+1)..(P+Q+R),1]  
		beta = J(cols(data),P+Q+R,.)		
		for(i=1;i<=cols(data);i++){
			X    = data[1..P,i]
			Y    = data[(P+1)..(P+Q),i]
			B    = data[(P+Q+1)..(P+Q+R),i]   		
			gx   = gvec[1..P,i]
			gy   = gvec[(P+1)..(P+Q),i]
			gb   = gvec[(P+Q+1)..(P+Q+R),i]  
			c=(X',-Y',B')
			lowerbd=J(1,length(c),0)
			upperbd=J(1,length(c),.)
			A1 = (-Xref',Yref',-Bref')
			b1 = J(cols(Xref),1,0)
			A2 = (diag(gx),J(P,Q,0),J(P,R,0))			
			b2 =-wx
			A3 = (J(Q,P,0),-diag(gy),J(Q,R,0))
			b3 =-wy		
			A4 = (J(R,P,0),J(R,Q,0),diag(gb))
			b4 =-wb                   
			Aie=A1 \ A2 \A3 \A4
			bie=b1 \ b2 \b3 \ b4
			lp.setCoefficients(c)
			lp.setInequality(Aie, bie)
			lp.setBounds(lowerbd, upperbd)
			theta=lp.optimize()
			if(lp.converged()==1){
				beta[i,.]=lp.parameters()			
			}
		}
       st_view(te=.,.,strname,flag)            
       te[.,.]=beta                       
   }  
   end                                                                    

                            
/////////

capture program drop nddfdual
program define nddfdual
    version 16
       gettoken invars 0:0, parse("=")
       gettoken var 0:0, parse("=")
       gettoken gopvars 0:0, parse(":")
       gettoken var 0:0, parse(":")
       syntax varlist [if] [in],gen(string) [wmat(name) gx(varlist) gy(varlist) gb(varlist) rflag(varname)]   
       local bopvars `varlist'
       local gv `gx' `gy' `gb' 
        marksample touse 
		markout `touse' `invars' `gopvars' `bopvars'
		if `"`rflag'"'==""{
            tempvar rflag
            qui gen byte `rflag' = 1
        }
		tempvar touse2
		mark `touse2' if `rflag'
		markout `touse2' `invars' `gopvars' `bopvars'	   
        local data `invars' `gopvars' `bopvars'
        foreach v in  `invars' `gopvars' `bopvars'{ //注意次序
            qui gen double `gen'_`v' = .
            local gens `gens' `gen'_`v'
        }
        local k: word count `invars'
        local l: word count `gopvars'
		local r: word count `bopvars'
	    if `"`gx'"'==""{      
	        foreach v in `invars'{      
	            tempvar g`v'      
	            qui gen `g`v''=-`v'      
	            local gx `gx' `g`v''      
	        }             
	    }      
	      
	    if `"`gy'"'==""{      
	        foreach v in `gopvars''{      
	            tempvar g`v'      
	            qui gen `g`v''=`v'      
            local gy `gy' `g`v''      
	        }             
	    }             
	    if `"`gb'"'==""{      
	        foreach v in `bopvars'{      
	            tempvar g`v'      
	            qui gen `g`v''=-`v'      
            local gb `gb' `g`v''      
        }             
	    }      
		
		local gv `gx' `gy' `gb'
		
	      
	   if "`wmat'"==""{    
           tempname wmat
		   matrix `wmat' =J(1,`k'+`l'+`r',1)
       }     
	
        mata: nddfdual("`data'","`touse'", "`touse2'",`k',`l',st_matrix("`wmat'")',"`gv'","`gens'")

end

cap mata mata drop nddfdual()  
   mata:   
       void function nddfdual( string scalar varname, 
                          string scalar flag,   
                          string scalar rflag,   
                          real   scalar P,  
                          real   scalar Q,  
                          real  colvector w,
                          string scalar gname,  //方向向量
                          string scalar strname )  
       {  
        data=st_data(.,varname,flag)        
        data=data'  
        dataref=st_data(.,varname,rflag)  
        dataref=dataref'  
        M=rows(data)
        R=M-P-Q
        Xref = dataref[1..P,.]
        Yref = dataref[(P+1)..(P+Q),.]
        Bref = dataref[(P+Q+1)..(P+Q+R),.]            
        gvec=st_data(.,gname,flag)  
        gvec=gvec'  
        class LinearProgram scalar lp  
        lp = LinearProgram()   
        lp.setMaxOrMin("min")  
		wx   = w[1..P,1]
		wy   = w[(P+1)..(P+Q),1]
		wb   = w[(P+Q+1)..(P+Q+R),1]  
		beta = J(cols(data),P+Q+R,.)		
		for(i=1;i<=cols(data);i++){
			X    = data[1..P,i]
			Y    = data[(P+1)..(P+Q),i]
			B    = data[(P+Q+1)..(P+Q+R),i]   		
			gx   = gvec[1..P,i]
			gy   = gvec[(P+1)..(P+Q),i]
			gb   = gvec[(P+Q+1)..(P+Q+R),i]  
			c=(X',-Y',B')
			lowerbd=J(1,length(c),0)
			upperbd=J(1,length(c),.)
			A1 = (-Xref',Yref',-Bref')
			b1 = J(cols(Xref),1,0)
			A2 = (diag(gx),J(P,Q,0),J(P,R,0))			
			b2 =-wx
			A3 = (J(Q,P,0),-diag(gy),J(Q,R,0))
			b3 =-wy		
			A4 = (J(R,P,0),J(R,Q,0),diag(gb))
			b4 =-wb                   
			Aie=A1 \ A2 \A3 \A4
			bie=b1 \ b2 \b3 \ b4
			lp.setCoefficients(c)
			lp.setInequality(Aie, bie)
			lp.setBounds(lowerbd, upperbd)
			theta=lp.optimize()
			if(lp.converged()==1){
				beta[i,.]=lp.parameters()			
			}
		}
       st_view(te=.,.,strname,flag)            
       te[.,.]=beta                       
   }  
   end  
                                           

	
***
use Ex5.dta,clear  
gen gl=-labor 
gen gk=-capital 
gen ge=-energy 
gen gy=gdp 
gen gb=-co2
matrix w=(1,1,1,1,1)
nddfdual labor capital energy = gdp:co2, gx(gl gk ge) gy(gy) gb(gb) wmat(w) gen(p_)
br

***	
	
	

*#第十五章 Stata编程的调试
*##15.1 mata编程的错误排查

	use Ex4.dta,clear  
	putmata data = (K L Y CO2),replace  
	mata:  
	        data = data'  
	        nx=3  
	        ny=1  
	        nb=1  
	        s=ny+nb  
	        Xref = data[1..nx,.]  
	        Yref = data[(nx+1)..(nx+ny),.]  
        Bref = data[(nx+ny+1)..(nx+ny+nb),.]  
	end 

/////////
	
	mata mata set matalnum on  
	use Ex4.dta,clear  
	putmata data = (K L Y CO2),replace  
	mata:  
	        nx=2  
	        ny=1  
	        nb=1  
	        s=ny+nb  
	        Xref = data[1..nx,.]  
	        Yref = data[(nx+1)..(nx+ny),.]  
	        Bref = data[(nx+ny+1)..(nx+ny+nb),.]  
	        q = LinearProgram()  
	        beta = J(cols(data),nx+ny+nb,.)  
        for(i=1;i<=cols(data);i++){  
	            X    = data[1..nx,i]  
	            Y    = data[(nx+1)..(nx+ny),i]  
            B    = data[(nx+ny+1)..(nx+ny+nb),i]      
            c=(Y',-X,-B')  
	            lowerbd=J(1,length(c),0)  
	            upperbd=J(1,length(c),.)  
	            A1 = (Yref',-Xref',-Bref')  
	            b1 = J(cols(Xref),1,0)  
	            A2 = (J(nx,ny,0),-nx*I(nx),J(nx,nb,0))  
            b2 =-1:/X  
	            A3 = (-s*diag(Y)+J(ny,1,Y'),J(ny,1,-X'),J(ny,1,-B'))  
	            b3 =-J(ny,1,1)  
	            A4 = (J(nb,1,Y'),J(ny,1,-X'),-s*diag(B)+J(ny,1,-B'))  
            b4 =-J(nb,1,1)                     
            Aie=A1 \ A2 \ A3 \ A4  
	            bie=b1 \ b2 \ b3 \ b4  
	            q.setCoefficients(c)  
	            q.setInequality(Aie, bie)  
	            q.setBounds(lowerbd, upperbd)  
            theta=q.optimize()  
	            if(q.converged()==1){  
                beta[i,.]=q.parameters()  
	              }  
	            }  
	            beta  
  end  

/////////
  
	use Ex4.dta,clear  
	putmata data = (K L Y CO2),replace  
	mata:  
	        nx=2  
	        ny=1  
	        nb=1  
	        s=ny+nb  
	        Xref = data[1..nx,.]  
	        Yref = data[(nx+1)..(nx+ny),.]  
	        Bref = data[(nx+ny+1)..(nx+ny+nb),.]  
	        q = LinearProgram()  
	        beta = J(cols(data),nx+ny+nb,.)  
	        for(i=1;i<=cols(data);i++){  
	            X    = data[1..nx,i]  
	            Y    = data[(nx+1)..(nx+ny),i]  
	            B    = data[(nx+ny+1)..(nx+ny+nb),i]      
	            c=(Y',-X',-B')                                            //word中X未转置  
	            lowerbd=J(1,length(c),0)  
	            upperbd=J(1,length(c),.)  
	            A1 = (Yref',-Xref',-Bref')  
	            b1 = J(cols(Xref),1,0)  
	            A2 = (J(nx,ny,0),-nx*I(nx),J(nx,nb,0))  
	            b2 =-1:/X  
	            A3 = (-s*diag(Y)+J(ny,1,Y'),J(ny,1,-X'),J(ny,1,-B'))  
	            b3 =-J(ny,1,1)  
	            A4 = (J(nb,1,Y'),J(ny,1,-X'),-s*diag(B)+J(ny,1,-B'))  
	            b4 =-J(nb,1,1)                     
	            Aie=A1 \ A2 \ A3 \ A4  
				bie=b1 \ b2 \ b3 \ b4  
	            q.setCoefficients(c)  
	            q.setInequality(Aie, bie)  
	            q.setBounds(lowerbd, upperbd)  
	            theta=q.optimize()  
            if(q.converged()==1){  
	                beta[i,.]=q.parameters()  
	              }  
	            }  
	            beta  
	  end

/////////	  
	  
	use Ex4.dta,clear  
	putmata data = (K L Y CO2),replace  
	mata:  
	        nx=2  
	        ny=1  
	        nb=1  
	        s=ny+nb  
	        Xref = data[1..nx,.]  
	        Yref = data[(nx+1)..(nx+ny),.]  
	        Bref = data[(nx+ny+1)..(nx+ny+nb),.]  
	        q = LinearProgram()  
	        beta = J(cols(data),nx+ny+nb,.)  
	        for(i=1;i<=1;i++){ // 修改为i<=1，只看第一个规划问题的约束条件  
	            X    = data[1..nx,i]  
	            Y    = data[(nx+1)..(nx+ny),i]  
	            B    = data[(nx+ny+1)..(nx+ny+nb),i]      
	            c=(Y',-X',-B')  
	            lowerbd=J(1,length(c),0)  
	            upperbd=J(1,length(c),.)  
	            A1 = (Yref',-Xref',-Bref')  
	            b1 = J(cols(Xref),1,0)  
	            A2 = (J(nx,ny,0),-nx*I(nx),J(nx,nb,0))  
	            b2 =-1:/X  
	            A3 = (-s*diag(Y)+J(ny,1,Y'),J(ny,1,-X'),J(ny,1,-B'))  
	            b3 =-J(ny,1,1)  
	            A4 = (J(nb,1,Y'),J(ny,1,-X'),-s*diag(B)+J(ny,1,-B'))  
	            b4 =-J(nb,1,1)                     
	            Aie=A1 \ A2 \ A3 \ A4  
	            bie=b1 \ b2 \ b3 \ b4  
	            q.setCoefficients(c)  
				q.setInequality(Aie, bie)  
	            q.setBounds(lowerbd, upperbd)  
	            q.getInequality()   //列印线性规划问题的不等式约束  
	            theta=q.optimize()  
	            if(q.converged()==1){  
	                beta[i,.]=q.parameters()  
	              }  
	            }  
				beta
	  end   

///////// 
	  
	cap mata mata drop deacrs()  
	mata:  
	    void function deacrs(string scalar d,   
	                         string scalar touse,  
	                         string scalar rflag,  
	                         real   scalar k,  
	                         string scalar res)  
	    {  
	        data    = st_data(.,d,touse)  
	        dataref = st_data(.,d,rflag)  
	        data    = data'  
	        dataref = dataref'  
	        class LinearProgram scalar lp  
	        lp = LinearProgram()                  
	        N1=cols(data)                      
	        N2=cols(dataref)                   
	        theta=J(N1,1,.)                     
	        M=rows(data)        
	   for(i=1;i<=N1;i++){                   
	        X=data[1..k,i]                      
	        Y=data[(k+1)..MM,i]                  
	        Xref=dataref[1..k,.]                
	        Yref=dataref[(k+1)..M,.]            
	        c = (1, J(1,N2,0))                  
	        lowerbd =., J(1,N2,0)               
	        upperbd = J(1,N2+1,.)               
	        Aie = (J(k,1,0),Xref \ Y,-Yref)     
	        bie = X \ J(M-k,1,0)                  
	        lp.setCoefficients(c)                
	        lp.setInequality(Aie, bie)           
			lp.setBounds(lowerbd, upperbd)       
	        theta[i]=lp.optimize()               
	    }  
	    te=1:/theta  
	    st_view(te=.,.,res,touse)  
	    te[.,.]=1:/theta  
	}  
	end                                                                        //提示不全
	/*
	variable M set but not used 并没有提示
	*/

*##15.2 Stata ado命令编程的错误排查

	cap program drop dea  
	program define dea  
	    version 16  
	    gettoken var 0:0, p("=")  
	    local input `var'  
	    gettoken var 0:0, p("=")  
	    syntax varlist [if] [in], te(string) [rflag(varname) rts(string)]  
	    marksample touse  
	    local output `varlist'  
	    if "`rflag'"==""{  
	        tempvar rflag  
	        qui gen `rflag' =1  
	    }  
	    local k: word count `input'  
	    qui gen double `teeff'=.  
	    if "`rts'"=="" local rts crs  
	    mata: dea`rts'("`input' `output'","`touse'","`rflag'",`k',"`te'")  
	end  
	  
	cap mata mata drop deacrs()  
	mata:  
	    void function deacrs(string scalar d,   
	                         string scalar touse,  
	                         string scalar rflag,  
	                         real   scalar k,  
	                         string scalar res)  
    {  
	        data    = st_data(.,d,touse)  
	        dataref = st_data(.,d,rflag)  
	        data    = data'  
	        dataref = dataref'  
	        class LinearProgram scalar lp  
	        lp = LinearProgram()                  
	        N1=cols(data)                      
	        N2=cols(dataref)                   
	        theta=J(N1,1,.)                     
	        M=rows(data)        
	   for(i=1;i<=N1;i++){                   
	        X=data[1..k,i]                      
	        Y=data[(k+1)..MM,i]                  
	        Xref=dataref[1..k,.]                
	        Yref=dataref[(k+1)..M,.]            
	        c = (1, J(1,N2,0))                  
	        lowerbd =., J(1,N2,0)               
	        upperbd = J(1,N2+1,.)               
	        Aie = (J(k,1,0),Xref \ Y,-Yref)     
	        bie = X \ J(M-k,1,0)                  
	        lp.setCoefficients(c)                
	        lp.setInequality(Aie, bie)           
	        lp.setBounds(lowerbd, upperbd)       
	        theta[i]=lp.optimize()               
	    }  
	    te=1:/theta  
	    st_view(te=.,.,res,touse)  
	    te[.,.]=1:/theta  
}  
	end    

///////// 

	set trace on  
	use Ex3,clear  
	dea K L = Y, te(eff) 
